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ABSTRACT 

In this paper we investigate the use of the vector potential as a means of maintaining the diver- 
gence constraint in the numerical solution of the equations of Magnetohydrodynamics (MHD) 
using the Smoothed Particle Hydrodynamics (SPH) method. We derive a self-consistent for- 
mulation of the equations of motion using a variational principle that is constrained by the 
numerical formulation of both the induction equation and the curl operator used to obtain 
Ph; the magnetic field, which guarantees exact and simultaneous conservation of momentum, en- 

Q . ergy and entropy in the numerical scheme. This leads to a novel formulation of the MHD 

force term, unique to the vector potential, which differs from previous formulations. We also 
demonstrate how dissipative terms can be correctly formulated for the vector potential such 
that the contribution to the entropy is positive definite and the total energy is conserved. 

On a standard suite of numerical tests in one, two and three dimensions we find firstly 
that the consistent formulation of the vector potential equations is unstable to the well-known 
SPH tensile instability, even more so than in the standard Smoothed Particle Magnetohydro- 
| dynamics (SPMHD) formulation where the magnetic field is evolved directly. Furthermore we 

■^J- . find that, whilst a hybrid approach based on the vector potential evolution equation coupled 

C^l ' with a standard force term gives good results for one and two dimensional problems (where 

. dA z /dt = 0), such an approach suffers from numerical instability in three dimensions related 

' to the unconstrained evolution of vector potential components. We conclude that use of the 

ON . vector potential is not a viable approach for Smoothed Particle Magnetohydrodynamics. 
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1 INTRODUCTION 

Magnetic fields are important at some level in almost every area of astrophysics and their effects are commonly accounted for in numerical 
simulations by solving the equations of magnetohydrodynamics (MHD). The solution of the MHD equations throws up a number of chal- 
lenges for any numerical method, most notably because of the divergence-free ("no monopoles") constraint on the magnetic field, which does 
not appear explicitly in the MHD equations but rather as an initial condition which, if satisfied initially, should remain satisfied for all time 



not appear explicitly in the MHD equa 
l T6thl200(]| : Price & Monaghan 2005). 



Attempts to solve the equations of MHD using the Smoothed Particle Hydrodynamics (SPH) method (for r ecent reviews see |Pricel2004 



Monagharl 2005) have a long and so mewhat tortured history, beginn ing with one of the earliest SPH papers bv lGingold & Monaghari" Il977l) 



though not seriously developed until Phillips & Monaghan! i 1985 ). The latter authors discovered that the equations of Smoothed Particle 



Magnetohydrodynamics (SPMHD) contained a catastrophic numerical instability (now known as the "tensile instability", Monaghan! 2000) 



when written in a form that con served momentum e xactly, albeit o nly in a certain regime (magnetic pressure greater than gas pressure). 



Despite detailed i nvestigation by Morris and Meglicki ( 1995), these problems meant that, apart from a few isolated applications (e.g 



Dolagetal.lll999l) , SPMHD did not find widespread adopt 



ion. 



More recently, progress has been made on a number of fronts, i n particular with regards to formulating dissipative terms in order 
to handle MHD shocks (Price & Monaghan 2004a, hereaft er IPaper j) a nd in form ulating the SPMHD equations self-consistently from 
a variational principle (Price & Monaghan 2004b, hereafter IPaper Ilh . In IPaper III we have furthermore derived the equations of motion 
accounting for terms related to the smoothing length gradients which are necessary for exact simultaneous conservation of both energy and 
entropy. There is also a reasonable consensus on good approaches to removing the tensile instability in SPMHD, using either formulations 
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proposed by Morris i 1996b or by B0rve. Omang & TrulsenI i 200 lb which forgo the exact conservation of momentum slightly in order to 
attain numerical stability, the latter method simp ly by subtracting off the ter m in the force equation which is proportional to the non-zero 
divergence. Application of the Monaghanl 1 2000b fix for the tensile instability l Paper 1 ) turned out to be unsuccessful in compressible flows 
i Priceil2004h . 



In Price & Monaghan (2005) (hereafter iPaper HID , an extensive investigation was made into methods for addressing the divergence 
constraint in the context of SPMHD, using either projection methods or hyperbolic cleaning schemes. However, whilst such schemes could 
be made to give good results on test problems (which can also be successfully run without any form of divergence cleaning), it was found 
that they perf ormed poorly for "re al" applications such as that of star formation where length and time scale s can change by many orders 
of magnitude i Price & Bate 2007 ). For this reason our attention shifted to an earlier formulation adopted by Phillips & Monaghanl i 1985b 



(though later discarded due to poor accuracy when calculated with the spatially constant smoothing lengths an d simple gr adient operators 



simple gradient operators 

used by these authors) wher eby the magnetic field was f ormulated in terms of the so-called "Euler potentials" dSternll 197cL 1976b (referred 
to as "Clebsch variables" by Phillips & Monaghari l 1985b . cue and (3e, where the magnetic field is expressed as 



B = X7a E x V/3 B . 



(1) 



The advantages are twofold - the first is that the divergence constraint is satisfied by construction (that is, taking the divergence of (TJ gives 
zero). The second is that the induction equation in ideal MHD written in terms of the Euler potentials takes a particularly simple form, namely 



= 0, 



(2) 



doLE _ dp E 
dt ~ ' dt 

which corresponds physically to the advection of magnetic field lines by Lagrangian particles dSte rnlll966b . 

Using the Euler potentials formulation together with the dissipative terms and force equation proposed in Paper Hit has meant that the 



SPMHD algorithm has been successfully applie d to a number of "real- world" probl ems, including neutron star merg ers (Price & Rosswog 
20060, star formation ( Price & Bate 2007, 2008) and the dynamics of spiral galaxies (Dobbs & Price 2008; Kotarba et al. 2009). Furthermore 
this has led to the development of at least two "magnetic-SPH" codes ( Rosswog & Pricell 20071 : Dolag & Stasvszvni 2008 ). the latter adding 
MHD to the widely used GADGET code for cosmological simulations (Springel 2005). 



widely used GADGET code for cosmological : 
However there are also important limitations to the Euler potentials approach, namely that the magnetic helicity 

A B dV, 



(3) 



where A is the magnetic vector potential and B is the magnetic field, is constrained to be zero by a simple consequence of the fact that 
A = q_bV/3_b (equivalently A = —(3e^cxe by a change of gauge) is exactly perpendicular to B, implying that A ■ B = 0. In practise 
this means that firstly it is simply not possible to represent certain fields using Euler pote ntials, as the y would become double valued. For 
example, one can easily represent either a toriodal or a poloidal field using Euler potentials Isternl 19761) but not a combination of bott0. The 
corollary is also that such complex fields cannot be created during the simulation. A better way of understanding this limitation in practise is 
to recognise that equation l[2}, since there is no time evolution of the potentials on the particles, represents a mapping of the magnetic field 
on the initial particle configuration at t = to a new arrangement on the particle configuration at some later time t, and can change the 
geometry of the field only insofar as a one-to-one mapping from the initial to the final particle distribution exists (i.e., the field can only be 
followed for < 1 dynamical time). Thus important physical processes such as winding up of fields by differential rotation are largely missed 
by the Euler potentials formulation (see lBrand enburg 2 0091 for some examples) and this motivates us to consider a more general approach. 

In this fourth paper we examine the question of whether a formulation of SPMHD based on the magnetic vector potential A can resolve 
the difficulties associated with the Euler potentials formulation whilst at the same time maintaining the divergence constraint on the magnetic 
field. In the vector potential formulation the magnetic field is given by 



B = Bi„i + Hext = V X A + Beart, 



(4) 



where 'Bint = V x A is the magnetic field due to internal currents i.e., those within the computation domain, and here we assume that Ji ex t 
is a time-independent externally applied field. The time evolution of A can be derived from the induction equation for the magnetic field 



^=Vx(vxB)-Vx (r?J), 



(5) 



where rj is the magnetic resistivity (77 = l/(u/io) where a is the conductivity and fj,o is the permeability of free space) and J = V x B//io 
is the current density. 'Uncurling' this equation, we have 



OA 

-=vxB-rjJ + V^, 



(6) 



1 Although it later turned out that there were problems in applying the Euler potentials formulation in the context of neutron star mergers. 

2 One possible way around this restriction is to note that one can co nstruct any given magne t ic fie ld by a linear combination of B's, where each B is 
determined by a separate set of Euler potentials, as done for example by Yahalo m & Lvnden- Bell ( 2008). 
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where <j> is an arbitrary scalar representing the freedom to choose a gauge. In Lagrangian form, and expanding B in terms of internal 
Bi nt = V x A and external H ex t (i.e., produced by currents outside the simulation domain) components this is given by 

dA 

_ = vxVxA + (v V) A + v x B? xt - r)3 + V</>. (7) 

Part of the difficulty in assessing the usefulness or otherwise of the vector potential is that differences between the Euler potentials 
and the vector potential onl y occur in three dimensions, since for one and two dimensional problems the two formulations are identical 
(e.g. lRosswog & Pricell2007r) . since in two dimensions the vector potential has only one component A z = oe which evolves according to 
dAz/dt — (for the Euler potentials V/3 = z in 2D). However a rigourous formulation of the SPMHD equations of motion has not been 
derived for either the vector or Euler potentials which has relevance to both formulations even in one and t wo dime nsions. 

The approach we take is to use a Lagrangian variational principle (similar to the approach taken in lPaper I] for the standard case) in 
order to derive the equations of motion in a manner that is constrained by the exact numerical formulation of both (0 and the numerical 
representation of the induction equation for the vector potential ( i]2.1b . This means that exact conservation of all physical quantities (linear and 
angular momentum and energy) is guaranteed in the resultant numerical equations provided that the appropriate symmetries (i.e., invariance 
to translations, rotations and time, respectively) are present in the Lagrangian and the equations used to constrain it. We find that this very 
powerful approach leads to a novel formulation of the force term which is already different to previous SPMHD formulations of the MHD 
force in one and two dimensions and indeed conserves momentum and energy exactly. We demonstrate that these symmetries are also 
respected in three dimensions provided an appropriate gauge choice is made in the dA/dt equation in order that it is Galilean invariant. 

Secondly, in Sj3]we show how dissipative terms should be constructed for vector potential SPMHD in order that total energy is con- 
served and that the second law of thermodynamics is obeyed, i.e., a positive definite contribution to the entropy results. These terms, which 
are derived independently of the equations of motion, differ from previous formulations of dissipative terms that have been used for the 
vector/Euler potentials in SPMHD. 

Finally, we examine the new vector potential force formulation and the dissipative terms on the suite of one and two dimensional 
test problems (Sj4j. Whilst the hope was that by constructing the SPMHD equations such that the divergence-free constraint was inbuilt, 
instabilities would not appear in the equations. However it turns out that the consistent formulation of the vector potential force has similar 
- in fact, much worse - problems with the tensile instability than even the standard conservative SPMHD force. Whilst we have managed to 
obtain reasonable results on a range of numerical tests with the consistent vector potential equations of motion, we find that a bet ter approach 



is to u se the vector potential in conjunction with a stable but non-conservative force such as those employed in Paper III and B0r ve et al 
(2001). The main practical improvement in this paper is therefore in the formulation of the dissipative terms. 



2 A CONSISTENT FORMULATION OF SPMHD USING THE VECTOR POTENTIAL 



2.1 Variational Principle 



We start from the Lagrangian for MHD, which in continuum form is given by (e.g. lNewcombll 19621) 



L = 



1 2 

2 PV 



pu 



1 

2(i 



B dV, 



(8) 



which is simply the kinetic minus the thermal and magnetic energies. The SPH Lagrangian is obtained, following Paper III and lMonaghan & Price 
(2001) by replacing the integral by a summation and the mass element pdV by the mass per SPH particle m, giving 



N 

~>sph — / TTlb 
6=1 



1 2 

2 Vb 



Ub(pi 



, 1 B{ 

< S b) - ~ 

2/io Pb 



(9) 



where v = x is the velocity, p is the density, u is the thermal energy per unit mass (in general a function of both density p and entropy s) 
and B is the magnetic field. 

The equations of motion can be derived using the Euler-Lagrange equati ons provided that all variables appearing in the La grangian 
can be expressed as a function of the particle coordinates and velocities (e.g. iMonaghan & Pricdl200 ll ; Price & Monaghanl 120070 . Whilst 
for hydrodynamics the density can be written directly as a function of the particle coordinates via the SPH density sum (which is an exact 
solution of the continuity equation), for MHD the magnetic field B (in this case the vector potential A) can only be written a s a function of 
the change in particle coordinates (i.e., we do not have an exact solution to the induction equation). In this case (as in lPaper IS) we can derive 
the equations of motion by perturbing the Lagrangian and specifying that the change in action is zero, i.e., 



8S - 



SLdt = 0. 



(10) 



where the variation SL is with respect to a small change in the particle coordinates <5x. Importantly conservation properties (e.g. momentum 
conservation) in this case will only follow provided that the respective symmetries (e.g. invariance to translation) are preserved in the 
numerical representation of the perturbation. 
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Perturbing the Lagrangian (Equation[9} with respect to a change in the position of particle a, i.e., <5x a , we have 



SL — m„ 



in b 



du b 



dp b 



Spb - 



2p 



Spb + 



Po Pt 



^ ■ 5 (p b B b ) 



(11) 



where we have expressed the perturbation for the magnetic field in terms of 8(pB) for reasons that will become clear. The equations of 
motion are obtained by using ( Hit in dlOt and integrating the velocity term by parts with respect to time, i.e., 

dv a 



) dt = [m a v a 



5x, 



dt 



5x a dt, 



giving (in tensor notation to avoid confusion of indices) 



' dt 



E 



nib 



Pb 8p b 

pI Sxi 



3 
2^0 



Spb 
Sxi 



+ 



1 H 5 (pbBl) 



Po p\ 



Sxi 



Sxi dt = 0, 



(12) 



(13) 



where we have used the first law of thermodynamics to write du/dp\ s = P/p 2 (see lPaper IIP . 

What remains is to express, as SPH summations over neighbou rs, the perturbations Spb and S(pbTib) taken with respect to <5x a . The 
derivation here is more complicated than that presented by Paper II because in this case we must express <5B in terms of SA (via the SPH 
expression of Equation f4]( and in turn, SA in terms of <5x (via the SPH version of the induction equation for the vector potential). We 
thus formulate the SPH expression of each of these equations in ij2.2| below, with the corresponding perturbations presented in £|2.4. 11 The 
equations of motion are derived in ^2.4.41 



2.2 SPH formulation 



2.2.1 Density sum 

We base our SPMHD for mulation for the vector po t ential on the variable smoothing length formulation of SPH presented by IPaper Tl ; 
Price & Monaghar] 1 2007 ) (see also Monaghar] 2002 . 2005 ). The density p is calculated on particle a from neighbouring particles via the 
summation 



p a = ^2m b W(\r a - r b \,h a ), 



(14) 



where W is the SPH kernel function, details of which are given in Appendix iBl Key to the variable smoothing length formulation is whilst 
the kernel function depends on the smoothing length h, h itself is defined as a function of the particle positions, expressed most conveniently 
as a function of the density sum itself, via the relation 



h a = Tj 

' Pa 

with derivatives 

dh a _ h a 
dp a vp a ' 



d 2 h a 



h a 

n2 



u+1 



(15) 



(16) 



dpi Pa 

where v is the number of spatial dimensions and r\ is a dimensionless constant specifying the smoothing length in terms of the mean inter- 
particle spacing. Equation J 1 5b in turn determines the "number of neighbours" in the SPH calculation. Unless otherwise indicated we use 
r\ = 1.2 throughout this paper, corresponding to ~ 58 neighbours in 3D for kernels with a compact support radius of 2h. 

The d ensity summation is theref ore a non-linear equation for both p and h that we solve using a Newton-Raphson scheme as described 
in detail in lPrice & Monaghanl 1 2007 ). Enforcement of the relationship between h and p is a necessary requirement for energy conservation, 
since hereafter we will assume that the smoothing length is differentiable with respect to particle position. 



2.2.2 V x A 

In principle we have a number of choices for the numerical formulation of both equation © and the induction equation. In practise these 
choices are constrained by the requirement that symmetries in the Lagrangian are preserved by the constraint equations. For example, in 
calculating 10 we can in principle use any of the SPH curl operators as discussed e.g. in Price ( 2004). The basic operation for a curl in SPH 
is given by 

(VX A) a = -Y—A h xV a Wab, (17) 
b P b 

where W a b = W^(|r a — r B \, h) is the SPH kernel. In principle the kernel used for the curl does not have to be the same kernel as used in the 
density summation, though in this paper we assume that this is the case. Since the time evolution of (and thus the perturbation to) A is not in 
itself invariant to translations in v (for example in the case of an external field, see below), we require a curl operator for V x A such that B 
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is invariant to the addition of an arbitrary constant to A such that the Lagrangian is Galilean invariant and therefore that the resultant force 
will conserve momentum. In the variable smoothing length formulation of SPH, this can be achieved using 



B a = (V X A) a - 



Em b (A a - A 6 ) x X7 a W ab (h a ) + B e 



where Q a is a normalisation factor related to the smoothing length gradient, given by 
Ma = 1 



dh a \ - dW ab 
— - > m b 



dh a 



(18) 



(19) 



which can be calculated alongside d!8t and the density sum d 14b . Note that we cou ld equally have chosen the alternative SPH curl formulation 
which has a 1 / p b inside the summation instead of the above (see e.g. Price! 2004) — equation dl8l ) has the advantage that it does not require 
prior computation of the density (which involves a summation over the particles) and depends only on the particle's own smoothing length 
(i.e., h a ). The symmetric formulation of the curl, i.e., 



(V x A) a 



-pa 22 m b 

b 



A a x VaW ab (h a )_ A b x V a W ab (h b ) 



Slapl M b pl 

is ruled out by the requirement that B be invariant to the addition of an arbitrary constant to A, as discussed above. 



(20) 



2.2.3 Time evolution of A 

For the induction equation we have freedom to choose not only the SPH formulation of the derivatives in Equation but also to choose 
an appropriate Gauge for the evolution of the vector potential (i.e., a choice of <f> in Equation0. Again here, our choice is constrained by 
physical requirements. Most importantly, in three dimensions in order to obtain momentum conservation in the equations of motion, we 
require that the induction equation is invariant to translations (i.e., is Galilean invariant). A gauge choice which achieves this, first suggested 
to us by Axel Brandenburg (2007, private communication) is to choose <f> = v • A, which leads to the induction equation in the form 

dA 

— = -A x (V x v) - (A • V)v + v x B ext - 77J. (21) 

It turns out that there are other good reasons for this choice of gauge - most notably that this in fact repre sents the correct low-speed (v « c) 
and magnetically dominated (E « cB) limit for electromagnetism 1 de Montignv & Rousseauxll2007 ). Written in tensor notation, equation 



can be expressed by 

dAi dv 3 j k 

~dT = ~ dx^ ijk v ext ~ r, Ji, (22) 

where etjk is the Levi-Civita permutation tensor and repeated indices imply a summation. Compare this with a "naive" gauge choice V0 = 0, 
which gives (from equation 

dAi j dAj j k 

-dT= V 3 Q^+^Be* t -vJi- (23) 

It is worth commenting that both and i2U are significantly more complicated than the evolution equations for the Euler potentials 
in 3D 0. The number of derivatives that require numerical evaluation for a Lagrangian code is also one more than would be required in an 
Eulerian scheme, since here we must compute not only the Eulerian term but also a "reverse advection" term to obtain the Lagrangian time 
derivative on the left hand side. Since advection terms are generally the source of the most difficulty in Eulerian codes - the main advantage 
of SPH is that these terms are not present - the accuracy with which these derivatives can be computed in practise on SPH particles is a 
concern. Whilst in the Galilean invariant gauge d21b the "reverse advection" term becomes a derivative of v rather than A, the same concerns 
apply. 

With regards to the SPH formulation of equation J22b we face a similar choice of SPH operators for computing both the curl and gradient 
terms as in ^2.2.21 In this case we are again constrained by the physical requirement that the SPH expression of i22l should be invariant to 
the addition of an arbitrary constant to v, which is achieved using a similar operator to that used for the curl of A. Neglecting dissipative 
terms (discussed separately in S0, we have 

dAi A] . ^ dWq b (ha) k 

~dT = na~p~a\ ^ ~ Vb) ^i— + (24) 

where in the above and throughout this paper, we adopt the convention that a, b, c, d refer to particle labels whilst i, j, k,l,m and n refer to 
vector/tensor components. Again, in principle we could also choose the form with l/p b inside the summation rather than the above. As in 
i]2.2.2l we choose the above expression because it does not require prior knowledge of the density to compute and can therefore be computed 
efficiently alongside the density summation. 
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2.3 Hybrid approach 



At this point a hybrid approach would be to compute the time evolution of the vector potential using equation J24b . calculate a magnetic field 
B using J 1 8b and then simply use this B in the equations of motion using (any of) the usual SPH expression(s) for the Lorentz force (e.g. 
Paper III : Paper IlJ) . for example the Morris 1 1996h formulation: 



dvl = 
dt 

where 

dxi 



E m " 

b 



Pa + ±B 2 a /fl dW ab (h a ) Pb + ^Bl/flodWabihb) 



dx l 



+ ■ 



dx' 



1 ^ (BiBj) b - {BiBj) a dW ab 
H mb — 



//n 



papb 



dxi 



dW ab (h a ) dW ab (h b ) 



dxi 



dxi 



(25) 



(26) 



Alternatively one can use the stable MHD force formulated by |B0rve et alj ( I2004I) . that is, where the source term |B(V ■ B) is subtracted 
from the conservative force, giving 



where 



Ml dW ab (h a ) 
n a pl dxi 



+ 



Ml 3 



dW ab {h b ) 



hpi 



dxi 



1 B x a \ ^ 

+ o — mh 

2 Mo V 



Bj dW ab (h a ) 

VaPi dxi 



+ 



dW ab (h b ) 



VbPl dxi 



1 

Mo 



B'B 3 



-Bo 



(27) 



(28) 



Indeed this is exactly the approach taken using the Euler potentials by e.g. Price & Bate i 2007 , 200 8h . The flaw in this methodology 
is that, since the equations of motion are not derived with the constraint of the numerical for mulation of the ind uction equation, there is no 
guarantee that total energy will be conserveqj (and indeed, using either the Morris force or the lB0rveetaIlj2OO4h approach, total momentum 
conservation is not guaranteed either, though the errors are quite small even for shock-type problems, see Price 2004). In this context total 
energy conservation means that 



dE 
~dt 



d_ 

dt 



E 



nib 



la, , 1 Bt 

-v b +u b + 

2 2 fiopi, 



E 



m b v b 



dv b 
~~dT 



du b 
~~dT 



1 Bl dpb 
dt 



2 Pi 



+ 



B b 

popb 



dB b 

dt 



0, 



(29) 



6 \ r^rv, b 

where dB / dt is the time derivative of Q8) which in turn involves the time derivative of A and thus our induction equation d24t (we derive the 
expression for dB/dt in Appendix[A]l. What is required is that the dv/dt term in the above is consistent with the dB/dt term resulting from 
the vector potential evolution. Needless to say, guaranteeing the conservation of energy in a vector potential approach is thus a complicated 
business, and one which is best achieved by following a variational approach. 

Whilst the hybrid approach works reasonably well for the Euler potentials (where the time evolution is zero according to equation [2}, 
for the more complicated evolution of the vector potential (Equation [24jl, using the induction equation to derive and thus constrain the MHD 
force term is more important. Furthermore, as we show below, this leads to a novel formulation of the Lorentz force in SPH which has 
not previously been considered and which has a number of interesting properties. In fact it should be possible to derive the corresponding 
formulation for the Euler potentials also, however we defer this to a future work. We compare the hybrid approach described above to the 
consistent vector potential formulation described below in the numerical tests presented in Sj4] 



2.4 Variational formulation 

2.4.1 Perturbations 

In order to derive the equations of motion from d 1 0t- J lib it remains to express the perturbations 8p b and 5( ppBb) in terms of the change in 
particle positions 5x a . The change in density is obtained by a perturbation of the density summation, giving (Paper II) 



dpi, 



^- ^2 m c (8x b - 5x c ) ■ X7 b W bc (h b ) 



which, when taken with respect to particle a, gives 
~- ^2m c (8 b a — 5 ca )\7bWbc(hb), 



5x a fib 



(30) 



(31) 



3 An important aside with respect to Eulerian codes is due here. Whilst in principle it is possible to enforce total energy conservation in any SPH scheme 
by simply evolving the total energy equation instead of the internal energy equation, if the system is not Hamiltonian all this does is push the errors to 
another quantity (for example the entropy in the case of hydrodynamics) and in practise would simply lead to negative pressures where 'total conservation' 
is violated. The very power of a Hamiltonian formulation of SPH is that exact and simultaneous conservation of all physical quantities is achieved (i.e., with 
zero dissipation) something which is never possible in a grid-based scheme. The caveat is that dissipation terms are then explicitly added to the SPH scheme in 
order to capture shocks and other discontinuities, the crudeness of which in practise often makes the scheme more dissipative than its grid-based counterpart. 
However there is no intrinsic dissipation in SPH. 
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where Sb a is the Kronecker delta referring to the particle labels. 

The derivation of the Lagrangian perturbation for the magnetic field from d!8b is a little more complicated and is given in Appendix lAl 
The result, including all terms relating to gradients in the smoothing length, is given by 



S(p b B b ) = — ^2 m c (A b - A c ) x [(5x 6 - <5x c ) ■ V] V 6 W bc {h b ) 



+ ^- E mc ( SAb ~ 5Ac ) x ^bW bc (h b ) + B ext 8p b 



Hi, 



B b 



n b 



e . Bfc ,intPb dh b ^ .. , . 

°Pb H k q — > m c [(dx 6 - dx c ) ■ VbJ 

it b ap b — 



dW bc (h b ) 
dhh 



(32) 



where we have assumed that B ex t is spatially constant (i.e., SB ex t = 0). The terms H and f , defined in Appendix lAl are higher order terms 
related to the gradient in the smoothing length which are necessary for strict conservation of energy - hence we retain them here - though 
they are generally expected to be negligible in practise. Taking the perturbation d32b with respect to Sx l a and using tensor notation gives 



Sxi 



tt b 



+ 



Sp b 



{5 ba — 5 ca ) 



d 

<-),! ' 



dW bc (h b ) 



dx{ 



-TT- > Die i 

fif> t—' V Sxi Sxi 



5A%\ dW bc (h b ) 
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The perturbation to the vector potential SA can be expressed as a function of <5x, from J24fr . by 



' 02; d J 77— Z 1 £kmnOX b I3 ext b . 

a xi 



(33) 



(34) 



At this point it is worth briefly pausing to examine the consequences of d32b for the equations of motion. In particular if we consider 
( 1321 ) and J34b together, ignoring the last two terms relating to smoothing length gradients, then one may observe that there are essentially 
three separate terms that will contribute to the force. We will refer to these as the '2D', 'external' and '3D' force terms. The 2D term arises 
from the first term in equation ( 132b and follows only from the SPH formulation of the curl used to construct B from A (i.e., Equation 1 18b. 
We refer to this as the 2D term because it is the only term which is present for purely two dimensional simulations where SA — (i.e., there 
is no external field). The 'external' term is present in the case of external fields and arises from the combination of the <5x x B ex ± term in 
d34b and the second term in Equation ( 132b . plus the third term from Equation ( 132b . Finally the 3D term arises from the combination of our 
choice both of gauge and SPH formulation of the induction equation (i.e., the first term of equation [34b and the curl operator via the second 
term in d32b . 

Whist the derivation of the 2D force is simply a matter of substituting the first term in d33b into d 1 3 b and simplifying, the external and 
3D force terms are more complicated since they involve first substituting the terms in d34b into the second term of d33b and in turn into d 1 3b . 

2.4.2 2.5D/B ext component 

For the 'external' force term, we can substitute the second term in d34b into the second term of d32b and add the third term from d32b to obtain, 
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Using the standard identity for the Levi-Civita tensor tjkitkmn = Sj n 5i m — 5j m Sin we have 
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(35) 



(36) 



The astute reader will note that Galilean invariance in the perturbation (and thus momentum conservation in the equations of motion) only 
follows if we make the simplifying assumption that the external magnetic field is spatially constant (and therefore independent of the particle 
positions). This is in fact physical since an external field with spatial gradients can impart momentum to the fluid and the equations in that 
case would not be expected to show Galilean invariance. In this paper we will deal only with spatially constant external fields, for which the 
above simplifies to (in vector form, where we have also substituted for 8p b using|30b 

S(p b B b ) ext = ext S~] m = (^ x i> _ ^ x c) ■ VbWbc(hb) - tj- E mc _ 5xc ) Bea:t ' ^bW bc (h b ). (37) 

\ L b \l b 

C C 

Taking the above perturbation with respect to 5x l a , the resultant term in d33b is given by 
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(38) 
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2.4.3 3D component 

The 3D force contribution is given by substituting the first term in d34b . taken with respect to Sx\, into the second term of J33b . giving 
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(39) 



The fact that the <5x is so deeply nested in the perturbation of B in the 3D case (i.e., via a summation for the curl of A [equation [32), and 
via a second summation for SA [ equation 1341 ), as we will see below, leads to a force term which is somewhat complicated to calculate. 
Nevertheless, it is a force term which preserves the basic symmetries that we asked for, namely momentum and energy conservation in the 
SPMHD equations. For example, using the naive or 'standard' gauge choice (Equation|23l> involves one fewer summations since the 5x is not 
nested inside a derivative in the perturbation to A. However, the perturbation is not Galilean invariant and can be straightforwardly shown to 
lead to a force that does not conserve momentum. 



2.4.4 Equations of motion 

Putting the perturbations J3 1 1 > and l |33l > [the second term of which has been expanded into d38b and I 139H into M2>\ we have 
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where we have collected the isotropic terms relating to smoothing length gradients into a single term by defining 
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(41) 



where H 1 and £ are defined in Appendix lAl Since the perturbation Sx 1 is arbitrary, upon simplification j40\ implies that the principle of least 
action is satisfied by the equations of motion in the form 
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where the current J k is defined according to 
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(42) 



(43) 
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noting that the swapping of indices in the permutation tensor ejki — > £kji leads to a change of sign in the corresponding term in equation 
l !42b . Although the curl operator in the above definition J43b is determined entirely by the variational principle, remarkably this is simply the 
standard SPH symmetric curl operator in the presence of a variable smoothing length (i.e., Equation 120b . The equations of motion can be 
expressed in a more compact notation by writing the isotropic, 2.5D/B e:E t and 3D terms in terms of a stress tensor and the 2D and 2D Vh 
terms in terms of differential operators, giving 
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,(44) 

(45) 
(46) 



P0 Slapa dp a ' 

Note that because the 2D terms cannot be represented by a stress tensor, S tJ does not represent the usual MHD stress tensor, since the Lorentz 
force in this case is composed of the divergence of the stress tensor plus the 2D terms. 

For the case of a constant smoothing length, the equations of motion simplify to (in vector notation) 
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(47) 



(48) 



At this point it is worth stepping back to consider the SPH formulations encapsulated by the force terms in ( 142b or equivalently, ( 144b 
or J47b . The most fundamental question is whether or not the magnetic force terms in the equations of motion derived above indeed are a 
representation of the Lorentz force when translated to the continuum limit. Since the proof is somewhat involved, the details and a translation 
of each of the terms into continuum form are given in Appendix lei Given that the equations of motion are indeed correct in the continuum 
limit, the following comments can be made about their numerical representation: 



(i) The isotropic term in Q42) is similar in form to the hydrodynamic SPH force and the usual isotropic MHD force in SPH. However, 
in this case the magnetic term is subtracted from the hydrodynamic pressure which implies that this term may be unstable to the clumping 
(tensile) instability caused by negative pressures in the regime where -^j^B 2 > P. 

(ii) The 2D and 3D terms present a novel formulation for the anisotropic magnetic force in SPH (strictly these terms also contain part of 
the isotropic force term - see Appendix Icl. The 2D term vanishes for constant A and is perpendicular to A yet remarkably both the 2D and 
3D terms conserve linear momentum exactly since they are antisymmetric in the particle index - implying that rn a dv a /dt = 0. 

(iii) Calculation of the 2D term involves use of the second derivative of the SPH kernel which is problematic using the cubic spline because 
the second derivative has discontinuous gradients. However this can be resolved using smoother kernels. 

(iv) For a purely external magnetic field the stress tensor b ecomes iS IJ = —PS 1 -' + l/po(Bl xt B J ext — S tJ \B 2 xt ) which is identical to 
the usual conservative SPMHD force term (as derived, e.g., bv lPaper Ilh . This part of the force, whilst conservat ive, is unstable to the tensile 
instability when the external magnetic pressure exceeds the gas pressure. The solution proposed by Paper III for this case was to simply 



subtract the constant term Bl xt Bl xt / po from the stress when a constant external magnetic field is imposed. 

(v) Calculation of the 3D term involves a triple summation over the particles — first, to calculate the density and simultaneously the 
magnetic field according to l |18b ; second to calculate J via ( 143 b or ( |48b ; and third, to calculate the force term. Use of this approach is 
therefore 1/3 more expensive than a standard SPMHD scheme. 



2.4.5 Instability correction 

As implied by items i) and iv) in the above discussion, we find in £|4]that the consistent variational formulation of the vector potential 
equations of motion ( 142b are in practise highly unstable to the tensile instability known to plague conservative SPMHD formu lations in the 
standard case where the induction equation for the magnetic field B is evolved l Phillips & Monaghanll 19851 ; Paper 1 ; Paper III ). Worse still, 
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we find that the equations are unstable for lower values of the magnetic field strength compared to the gas pressure (i.e., higher plasma /3) 
than in the standard case, consistent with our conjecture in item i), above. That is, we find that the 2D term in J42l > does not provide any 
significant stabilising influence over the negative stress arising from the isotropic term. 

For this reason we consider below ways of correcting the force term such that the net stress is positive in order to stabilise the formulation 
against the tensile instability. An obvious approach is to simply revert to the hybrid scheme discussed in i )2.3l Indeed in the end (SJ5J we 
conclude that this is the best approach to implementing the vector potential in SPH, though one is left with the same problems that the 
consistent formulation given in i]2.4l was constructed to solve, namely that energy is not conserved exactly and that there is no (Hamiltonian) 
constraint on the overall evolution of the system (see discussion in ^2.3b . To this argument it may be countered that the consistent approach 
is no better in this respect once the correction terms below are added. 



In this paper, we stabilise the vector potential formulation we use the method proposed by |B0rve et alJ fcOOl"). namely adding the "source 
term" — B(V ■ B)/p to the (conservative) force, giving 



—r- = — ; B > rrib 



Bi dW ab (h a ) B\ dW ab (h b ) 



dxi 
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(49) 



This method violates both the conservation of momentum and energy, but only to the extent that the numerical estimate of V ■ B according 
to the above is no n-zero. Whilst the proof that this correction term indeed stabilises the standard SPMHD equations has been given by 
B0rve et alj J2004h . it is unclear why it should also work for the vector potential formulation. Empirically, we find that subtracting the 
B(V ■ B) term can indeed stabilise the vector potential force, but only in a limited range of circumstances, the limitations of which are 
unclear. For example, the ID shock tube problems ( £|4.1t are stabilised effectively by this method and similarly the circularly polarised 
Alfven wave in 2D ( ij4.2t . However the 2D Orszag-Tang Vortex ( £j4.3fc remains unstable even with the correction term added. Ideally a full 
stability analysis of the vector potential equations of motion J42t should be carried out, though such a task is well beyond the scope of this 
paper. 



A more severe alternative would be to use the original method of Philli ps & Monaghanl l l 19851) . namely to correct the stress by subtracting 
the maximum value over all the particles, 



Smax 



(50) 



Whilst this method conserves momentum but not energy, the correction to the stress can become arbitrarily large. For the Orszag-Tang Vortex 
problem ( i)4.3b we find that the stress correction required to stabilise the solution starts to produce unphysical features in the solution and is 
therefore an unacceptable alternative. 



3 DISSIPATIVE TERMS 

In Papers I and III the need to introduce dissipative terms in the magnetic field in order to account for discontinuities in B was discussed. For 
the magnetic field this means adding an artificial resistivity term. The key constraints on deriving such a term are that it should i) conserve 
total energy and ii) result in a positive definite increase in entropy - or equivalently - thermal energy. [Paper J used these constraints to derive 
appropriate artificial resistivity terms for the standard SPMHD approach (i.e., using B orB/p in the induction equation). 



3.1 Resistivity using the vector potential 



The formulation of resistivity in the vector potential formulation is considerably simplified since the derivatives of r\ do not enter the evolution 
equation for the vector potential (c.f. equation |2"0. The appropriate dissipative term in the vector potential evolution is therefore given by 



dA a 

dt 



(51) 



where the resistivity r\ a and an SPH expression for J remain to be defined. 

The constraint of total energy conservation is expressed by equation J29t . Since the magnetic dissipation (physically) does not enter the 
equations of motion nor the continuity equation, we are left with the requirement that 

' du a \ B a / dB a ' 



dt 



+ ■ 

diss MOPa 



dt 



= 0. 



(52) 



The reader should note that — whilst we do not consider dissipative terms as part of the derivation of the equations of motion from the 
Lagrangian in i)2.4l — the above would be equivalent to stating that the contribution to the Lagrangian from the perturbation to the magnetic 
evolution is exactly balanced by the contribution from the perturbation to internal energy (from an increase in entropy) in d I lb - thus having 
no effect on the equations of motion. Writing dB/dt in terms of dA/dt using JAlb gives 
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where the above expression is determined by our choice of SPH operator for V x A (see Appendix [At and we have dropped the subscript 
diss from the time derivatives assuming we are referring to the magnetic dissipation terms only. Using J5U in the above gives 



Ed.Ua 



^ m b (rjaJa - r/ b J b ) X X7 a W ab (h a ). 



(54) 



Adding half of this term to half of the same term with the indices a and b exchanged, using the antisymmetry of the kernel V bWba(hb) 
-V ' a Wab(hb) and the vector identity B ■ (r/J x VW) = -??J(B x X7W) gives 
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Swapping summation indices in the second term and combining the two terms, we have 
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giving the contribution to the thermal energy equation for particle a as 

B a „. , Bb 
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(55) 



(56) 
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(58) 



If we define J in d51t using the symmetric curl (eauationl20b. then we have simply 

du a \ _ r] a Jq 

dt ) drss Pa 

which is exactly the continuum expression. Furthermore the dissipation is guaranteed to be positive definite so long as J in gil l (and 
is calculated using the symmetric curl operator. 

As an aside, it is interesting to note that, as in the equations of motion, in the dissipation we are required to use the symmetric curl 
operator for J to obtain energy conservation, where in this case our only choice of SPH formalism was to specify the operator used in 
B = V x A. The reason for the appearance of the s ymmetric curl is because the curl operators in dl8t and J20I ) form a conjugate pair. This 
conjugacy in SPH operators has been noted earlier by Cummins & Rudmar] i 1999b in the context of divergence cleaning, though not with the 
variable smoothing length formulation. 

The disadvantage of the symmetric curl is that it can give a poor representation of J if the particles are disordered. In particular on the 
shock tube tests (fj4]l we find problems using the symmetric curl for the dissipation in combination with the consistent SPMHD equations of 
motion derived in ^2.4.41 though using exactly the same dissipation in combination with either {25} or {27} gives good results. A compromise 
approach is to use the usual curl operator (as in equation|4} for J in J51t and {58}. Provided that J 2 is used in {58} (rather than using|57b 
then the dissipation will be guaranteed to be positive definite, however energy conservation will only be approximate because the energy- 
conserving expression is given by {57}. Alternatively energy conservation can be enforced by using {57}, whereby positivity of the dissipation 
is only guaranteed so long as the J estimate used in {5J~} and the symmetric curl estimate of J given by the summation in ( 1571 ) have the same 
sign. 



3.2 Choosing 77 corresponding to an artificial resistivity 

The remaining issue is to formulate the resistivity parameter r\ appropriately for an artificial resistivity, that is where the resistivity acts only 
on the smallest scales in the calculation to diffuse discontinuities in the magnetic field. We propose the following: 

T} a = a B v a A h a , (59) 

where as is a dimensionless factor of order unity, v\ is the Alfven speed and h is the particle's smoothing length. An alternative which is 
second order in the smoothing length would be to use 



rj a = a B \ h a , (60) 

V 

where J is the current density. This gives a resistive diffusion that is nonlinear in the smoothing length and responds only to large gradients 
in B. In general we find ([{4} that {60} gives insufficient dissipation at discontinuities, though perhaps some combination of {59} and {60} 
could be a reasonable compromise. 
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3.3 Comparison with previous formulations 



The above formulation for the vector potential differs considerably from that required in the usual SPMHD approach dPaper J) since r\ 
cannot be defined using the pairwise signal velocity bet ween particles and t heir m utua l separation. The de rivation given above also shows 
that the naive dissipative terms previously formulated bv lRosswog & Pried (120070 and lPrice & B ate (2007) for the vector/Euler potentials, 
in volving the pairwise sign al velocity, have no guarantee of positive definite contributions to the thermal energy. The dissipation term used 
by Rosswop & Price! J2007 ) was given by 



dA a 

dt 



Ei ri b , a 

■^a B v si g (A a 

, Pab 



At) F ab 



(61) 



where F a b = \ [F a b(h a ) + F a b(hb)] is the scalar part of the kernel gradient term (see Appendix|B](. Following the analysis given above, it 
can be shown that the appropriate term to be added to the thermal energy equation in order to conserve total energy is given by 



du a \ 1 v ^ 

dt I - 



OtBV. 



m b — — - (Jo - Jfc) ■ (A a - A b )F ab 
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(62) 



The dot product between J ab and A ao has no guarantee of being positive definit e, and indeed in num erical tests on the sh ock tube problems 
discussed in Sj4]we find that nega tive therm al energies can result. Instead lRosswog & P rice (2007) and lPrice & Bate! ( 120071) add the term used 
in the usual SPH formulation bv lPaper III i.e., 



du a 

~dT 



m b r-^ (Bo 

POPab 



BbY F a 



(63) 



which is positive definite, but when combined with ( 16 1 1 has no guarantee of conserving energy. Thus, this approach should be discarded in 
favour of the correct formulation given by equations d51b and d58b . Note that J5 1 1 > and ( 1581 ) should be used regardless of whether the equations 
of motion derived in i)2.4,4| are implemented. That is, even using an alternative representation for the force terms such as d25b , the dissipation 
for the vector potential should still be implemented using equations d5U and d58b . the expressions for J in which are determined entirely by 
the choice of curl operator in B = V x A (see i]2,2,2b . There are similar implications for the dissipative terms used with the Euler potentials, 
though these will be discussed elsewhere. 



4 NUMERICAL TESTS 



We have implemented the vector potential formulation of SPMH D into the iV-dimensional SPH code (hereafter, NDSPMHD) that we have 
previously used to test the standard SPMHD formulation in Price ( 2004) and Papers I — III. The code evolves the SPMHD equations using a 



standard leapfrog predi ctor-corre ctor scheme, where the vector potential A is evolved alongside the velocity field. The timestep is controlled 
globally as described in lPaper IIIl where unless specified we use a Courant factor of C CO ur = 0.30. 



Use of the vector potential compared to the standard variable-smoothing length SPMHD scheme dPaper IIB) involves the following 
changes to the code: 

(i) During the iterated loop over the particles to calculate both density, smoothing length and Q, from (19), calculate B using ( 118b and the 
smoothing length gradient terms H and £ using dA4b and dA9b respectively. Note that H and £ can be combined with B at the end of the 
summation loop to construct £ according to d41b and thus stored only as a single scalar variable. 

(ii) In three dimensions an extra loop over the particles is required to calculate the current J using the symmetric curl (equation [43j. 

(iii) The force is calculated in the main loop according to ( 142b . alongside which the time derivative of the vector potential is evaluated 
using d24b . with dissipation according to d51b and d58b . 

(iv) Where external fields are used, boundary (ghost) particles require that the boundary value of the vector potential evolves according to 
the second term in d24b . 

In the following sections we examine the performance of the consistent formulation of the vector potential on a range of test problems 
that are commonly used to test (both grid based and SPH) MHD codes. The tes ts are identical to those de scribed in Papers I — III apart from 
the 3D version of the Orzsag-Tang Vortex which has also been considered by Dola g & Stasvszvr] J2008j>. Our first aim is to compare the 
energy-conserving or "consistent" formulation derived in £]2.4|bofh to the standard S PMHD scheme ( Paper III) but also to the more naive 
hybrid approach (i ]2.3b used in previous papers IPrice & BatdboollRosswog & PricdboOTt) . where in this paper we use d27b though results 
are similar with the Morris force. For this purpose, one and two dimensional problems suffice since the equations of motion d42b already 
differ from the hybrid approach using d27b or d25b in one dimension. Shock tube problems are particularly valuable for assessing the role 
of dissipative terms, for which the formulation used in this paper also differs from those used previously (see ^3. 3b . as well as the effect 
of the "2.5D" terms in d42b relating to an external magnetic field. The Orszag-Tang Vortex in 2D, though the most complicated in terms of 
dynamics, is the simplest test problem in terms of implementation, since the only non-zero terms present from d42b are the isotropic and 2D 
terms. 

Our second aim is to examine the accuracy of the vector potential in general, using either the energy-conserving ("consistent") or hybrid 
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Figure 1. Results of the Brio & Wu (1988) shock tube problem using the consistent vector potential formulation 142 \ corrected for stability using )49t . using 
the quintic kernel instead of the cubic spline and replacing the symmetric curl operator required in the MHD dissipation terms with a more accurate but 
n on-conservat ive estimate. Artificial viscosity, conductivity and resistivity have been applied, the first of these using a switch, the second using the formulation 
of Price (2008) and the third as in fusing a B = 0.75. 



approach. Testing of the vector potential evolution equation l !24t as distinct from the Euler Potentials (SjT} requires a three dimensional test 
problem, for which we consider a 3D version of the Orszag-Tang Vortex. Since we find that this problem cannot be effectively stabilised 
using the consistent formulation, we consider it only with the hybrid approach. 



4.1 1.5D shock tube problems 

The shock tube descri bed by Brio & Wul ( 19881) is perhaps the most widely used test problem for MHD codes. We set up the problem 



exactly as described in lPaper Il with no smoothing of the initial conditions and using approximately 800 equal mass particles in the domain 
x = [—0.5, 0.5]. Conditions to the left of the shock are given by (p, P, v x ,v y ,B y ) = [1, 1, 0, 0, 1] and to the right by (p, P, v x ,v y , B y ) — 
[0.125, 0.1, 0, 0, —1] with B x — 0.75 and 7 = 2.0. For the vector potential this means that we set Ji ex t = [B x , 0, 0] and A z = —B y x 
initially. 

The results using the consistent formulation of the vector potential force (equation [42} are shown in Figure Q] and may be compared to 
the exact solution given by the solid line. There are numerous issues. To produce a reasonable solution at all on this problem, we have had 
to: 

• Stabilise the force term using d49l >, 

• use the quintic kernel dB3t instead of the cubic spline dB2b , and 

• use a more accurate J estimate in the dissipation term d51t rather than the symmetric curl. 

The solution shown in Figure[T]also requires a relatively high dissipation parameter (q_b = 0.75) that is not applied using a switch. The 
fact that we have not used the symmetric curl means that the dissipation term does not conserve energy exactly — we have instead ensured 
that the contribution to the entropy is positive definite (see i)3.1b . We also checked whether or not simply subtracting the constant external 
component from the stress tensor as discussed in point iv) of i)2.4.4l would stabilise the result instead of {49} but found this not to be the case. 

The results using the hybrid approach (i.e., using ([27} for the force) are shown in Figuref2]and were obtained with far fewer tweaks — 
that is, using the cubic spline, ob =0.5 and the dissipation applied using the symmetric curl for J as required for energy conservation. The 
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Figure 2. Results of the lBrio & wl Jl988h shock tube probelm using a hybrid but non-conservative approach. In this case we have been able to use the usual 
cubic spline kernel and the correct curl formulation to ensure that the MHD dissipation is both positive definite and energy-conserving. Artificial viscosity, 
conductivity and resistivity have been applied as in Figure[T] though with a slightly lower resistivity parameter (ceg = 0.5). 



remaining issues with over-smoothing of the fast rarefaction are similarly present in the standard SPMHD formulation (see Paper II) and are 
mainly due to the fact that the resolution is very low in this region due to the density jump. The oscillations around the slow magnetosonic 
shock can be calmed further by increasing the resist ivity parameter, though a t the expense of further smoothing the fast rarefaction. The 
results in this case are similar to those presented by [Rosswog & Price! ( 120070 which is to be expected since the only difference is in the 
correct f ormulation of the dissipative terms used here. The results in Figures[T]and|2]are an improvement over tho se present ed in lPaper II and 
Paper III although this is mainly because of a better understanding of how to apply dissipation terms developed in Paper III rather than being 
an intrinsic i mprovement due to the use of the vector potential. 

As the lBrio & Wul dl988l) problem is more difficult in SPH because of the density cont rast, we consider tw o further shock tube problems 
in this paper, both of which have also been used to test the standard SPMHD formulation l Pricell2004 JPaper J) . 

The second shock tube illustrates the formation of seven discontinuities in the same problem (Figure[3j. The Riemann problem is set up 
with initial conditions to the left (x < 0) of the shock given by (p, P, v x , v y , v z , B y , B z ) = [1.08, 0.95, 1.2, 0.01, 0.5, 3.6/(4tt) 1/2 , 2/(4tt) 1/2 ] 
whilst to the right (x > 0) (p, P,v x ,v y ,v z , B y , B z ) = [1, 1, 0, 0, 0, 4/(4tt) 1/2 , 2/(4tt) 1/2 ] with B x = 2/(4tt) 1/2 and 7 = 5/3. Using 
the vector potential we set B Mt = [B x , 0, B z ] and A z — —B y x initially. Since the velocity in the x-direction is non-zero at the boundary, 
we continually inject particles into the left half of the domain with the appropriate left state properties. The resolution therefore varies from 
an initial 700 particles to 875 particles at t — 0.2. The results are shown in Figure[3]at time t = 0.2, using the consistent vector potential 
force formulation. As with the Brio-Wu problem, we find that it is necessary to correct the force term according to J49t in order to obtain 
stability, though in this case the solution is obtained satisfactorily using the cubic spline kernel, the symmetric curl in the dissipation and a 
low resistivity parameter as = 0.15. Figure[5]may be compared with Figur e 4 in lPaper IL though as in the Brio-Wu problem the differences 
compared to that paper are more due to the improvements made in lPaper ill rather than anything specific to the vector potential. Results with 
the hybrid formulation (i )2.3b are similar. 



The final shock tube problem we consider has initial conditions to the left of the shock given by (p, P, v x , v y , B y ) = [1, 20, 10, 0, 5/(47r) 1/ ' 2 ] 



and to the right by (p, P, v x , v y , u y 



B y ) = [1, 1, -10,0, 5/(4tt) 1/2 ] whhB x = 5.0/(4tt) 1/2 and 7 = 5/3. The vector potential is set up using 
B ea;t = [B x , 0, 0] and A z = —B y x initially. The results computed us ing the consistent vector potential formulation are shown at t — 0.08 
in Figure \4\ and may be compared with the exact solution taken from Dai & Woodward 1994 ) given by the solid line. We have used 800 
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Figure 3. Results of the adiabatic MHD shock tube problem showing the formation of seven distinct discontinuities related to the propagation of MHD waves. 
The solution has been computed using the consistent vector potential force, corrected for stability with {49}, using the standard cubic spline kernel and an 
artificial resistivity parameter ctg = 0.15. 



particles in the domain. Figure|4]may be compared with Figures 4.18 and 4.19 of iPricd ( 120041) for the standard SPMHD formulation. For 
this problem, apart from applying the stability correction J49t we have nonetheless used the cubic spline kernel, though with the consistent 
formulation satisfactory results could not be obtained using the symmetric curl in the dissipation and thus we have reverted to the more 
accurate (asymmetric ) J estima te, forgoing exact energy conservation. H ere we have a pplied both artificial viscosity and resistivity using the 
switches discussed in Paper III, and artificial conductivity as described in Price i 2008) with a u — 1. 



4.2 2.5D Circularly Polarised Alfven wave 

The circularly polarised Alfven wave is an exact, non-linear solution of the MHD equations. It is particularly useful as a test problem as it 
allows one to compute the evolution of a non-linear wave of arbitrary amplitude indefinitely, since the wave does not com press the ga s and 



therefore does not steepen into a shock. The parameters for the test problem used here are identical to those described by Tothl 1 2000l) for 



Eulerian codes and the setup for SPH is identical to that described in Paper III for the standard SPMHD scheme except that here we set up 
the magnetic field in terms of the v ector pote ntial. 



Whilst the reader is referred to 



> lPaper I II for further details, we briefly recap the setup parameters: The wave is setup in a two dimensional 
domain with a unit wavelength along the direction of propagation (ie. in this case along the line at an angle of 30° with respect to the x-axis). 
The initial conditions are p = 1, P = 0.1, v\\ = 0, B\\ = 1, v± = B± — 0.1 sin (27rr||) and v z = B z — 0.1 cos (27rr||) with 7 = 5/3 
(where r\\ = x cos 9 + y sin 9). The x— and y— components of the magnetic field are therefore given by B x = B\\ cos 9 — B± sin 9 and 
By = By sin# + B± cos6> (and similarly for the velocity and vector potential components). Conversely, Bu = B y sin 8 + B x cos 9 and 
B± = By cos 9 — B x sin 9. 

For the vector potential we set up two components, A± — 0.1[sin (27rrii)]/27r and A z = 0.1[cos (27rrii)]/27r, together with the Bi\ 
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Figure 4. Results using the consistent vector potential formulation of the MHD shock tube problem with a strong compression and weaker shearing disconti- 
nuities. Here we have added the correction to the force term for stability, the standard cubic spline kernel and an more accurate but non-conservative curl in 
the MHD dissipation terms which are applied using a switch. The results may be compared with the exact solution given by the solid line. 



component added as an external field. As the external field is the dominant component of the magnetic field, this is mainly a test of the 2D 
and 2.5D/B ex t terms in J42b together with the evolution of the magnetic field according to ( 124b . Furthermore, since the magnetic pressure 
exceeds the gas pressure, this test is unstable to the SPH tensile instability in the standard formulation (Paper I). Using our consistent vector 
potential formulation (equation 142b we expect the force to be unstable in the regime where P > 3/2B 2 (see note i. in i ]2.4.4t . In order 
to assess whether or not this was true in practice we computed a series of tests adjusting the value of the gas pressure (since the test is 
non-compressive and the wave travels at the Alfven speed, adjusting the gas pressure should have no effect on the results). Indeed we find 
that, in the absence of instability corrections the test is indeed unstable to the tensile instability and that this occurs for all simulations where 
P £ 3. 



4.3 2D Orszag-Tang vortex 

The compressibl e version of the Orszag & Tang J 19791) v ortex has also been widely used as a test problem for both grid based an d SPH 
MHD codes (e.g. lRvu et alll995l ; lpai & Woodwardll998 jLondrillo & Del ZannJ200d ; lT6fl]|200d : |p"aper IIltlRosswog & Pricell2007l) .The 
setup we use, identical to that in Irvu et alJ Jl995n and lLondrillo & Del Zannal (I2000h consists of an initially uniform density, periodic 
lxl box given an initial velocity perturbation v = vo[— sin (27rj/), sin (2irx)] where no = 1. The magnetic field is given a doubly 
periodic geometry B = Bo\— sin (2-iry), sin (47ra;)] where Bo = l/v47T> which for the vector potential is achieved by setting A z — 
—Bo/n [i cos (2ny) + j cos (47ra;)1 initially. No external fields are present for this problem. The flow has an initial average Mach number 
of unity, a ratio of magnetic to thermal pressure of j3 — 10/3 and ratio of specific heats 7 = 5/3. The initial gas state is therefore 
P = 5/3B 2 = 5/(12tt) and p = jP/v = 25/(36tt). In this paper we have used 128 2 particles, initially placed on a cubic lattice, which 
is sufficient to demonstrate energy conservation and the problems with numerical stability faced b y the energy-conserv i ng fo rmulation. 

We have considered the problem in its original 2D form and also in a 3D geometry following Dolag & Stasvszvnl 1 20081) . As discussed 
above the Orszag-Tang Vortex in 2D is the simplest test problem for the vector potential in terms of implementation, since it is purely two 
dimensional — in the magnetic field (i.e., no external fields) as well as the spatial coordinates. As there is no inflow at the boundaries it is 
also a good test for checking energy conservation in the code. 
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Figure 5. Results of the circularly polarized Alfven wave test at t = 5 (corresponding to 5 wave periods). The plots show the perpendicular component 
of the magnetic field vector B± = B y cos 9 — B x sin 9 for all of the particles, projected against a vector parallel to the direction of wave propagation 
rii = xcos9 + ys'm9 (where 9 = 30° in this case). The SPMHD results are shown at five different resolutions which are, from bottom to top, 8 X 16, 
16 X 32, 32 X 64, 64 X 128 and 128 X 256. The exact solution is given by the solid line. The left panel shows the results using the conservative vector potential 
formulation with the constant external stress subtracted whilst the right panel shows results using a 'hybrid approach' - that is, evolving the vector potential 
but with a standard SPMHD force term. Artificial viscosity was applied using a switch, though no artificial resistivity. The results are indistinguishable except 
that the highest resolution calculation in the left panel starts to show noise in the particle distribution at late times related to the excitation of compressible 
modes, most likely related to the poor accuracy in the kernel second derivatives. 



We find that the calculation using the energy-conserving vector potential formulation becomes unstable to the tensile instability at 
relatively early times (t > 0.1), manifested by the usual symptom of particles clumping along the field lines in the low density, high 
magnetic pressure regions. The development of the instability is clear by t = 0.2 as shown in the density field in the left panel of Figure[6] 
Worse still, adding the correction term l |49t does not stabilise the instability for this problem. In fact the only fail-safe method we have found 
of removing the instability in this case is to subtract the maximum value from the stress tensor as in d50t— and we find that doing so in this 
case produces unphysical artefacts in the solution. Thus we are unable to produce a stable and accurate solution to the Orszag-Tang Vortex 
using the consistent vector potential formulation. 

Using the hybrid approach the solution has no such difficulties (right panel), identical to the results presented in iRosswog & Price 



(2007). What is remarkable though is that despite the severe numerical instability present with the consistent formulation, energy is in fact 
much better conserved than using the hybrid scheme, as demonstrated in Figure[7J This figure shows the evolution of the total energy (left) 
and magnitude of the linear momentum (right) for the two calculations shown in Figure [6] and also for two calculations using a smaller 
Courant factor. For the consistent formulation, it can be seen from Figure[6]that energy is conserved exactly, i.e., to timestepping accuracy 
— meaning that the energy conservation can be improved to arbitrary precision by decreasing the Courant factor. By contrast, decreasing the 
timestep further in the hybrid case produces essentially no change in the energy conservation, indicating that the non-conservation derives 
from the SPH scheme itself. Indeed we find that this is a rigorous test of the implementation of the energy-conserving formulation in the 
code, since even by neglecting the high-order terms relating to the smoothing length gradients discussed in Appendix lAl we find that the total 
energy rises rapidly to unphysical values once the instability sets in. 

In fact, a very g ood solution to th e 2D O rszag-Tang Vortex problem can be obtained using the hybrid approach, as has already 
been demonstrated by Rosswog & Pried J2007). We do not feel that it is necessary to repeat those results here, referring the reader to 



Rosswog & Pried ( 120071) for details. Instead, in this paper we skip directly to the three dimensional version of the problem since we are 
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Figure 6. Density in the 2D Orszag-Tang vortex problem at t = 0.2 in a 128 2 particle calculation, showing development of the numerical tensile instability in 
the energy-conserving formulation of the vector potential (left), compared to a hybrid approach using a stable but non-conservative force formulation (right). 
The instability develops in the low density regions of the flow where the magnetic pressure exceeds gas pressure. Whilst the standard SPMHD formulation 
evolving B suffers from similar problems in the j3 < 1 regime when a conservative force is used, the onset of the instability for the conservative vector 
potential formulation occurs at lower magnetic field strengths (j3 < 3). 
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Figure 7. Energy conservation in the 2D Orzsag-Tang problem, comparing 128 2 particle calculations using the consistent formulation (solid/black, dashed/red 
lines) to the hybrid approach (dot-dashed/green, dotted/blue lines), for two different settings for the Courant factor in the timestepping (C CO ur as indicated). 
The corresponding density field at t = 0.2 is shown in Figure[6] Despite the strong numerical instabilities present in the solution with the consistent formulation 
(Figure[6), energy is conserved exactly to the accuracy of the timestepping scheme and linear momentum (right panel) is conserved to machine precision. Total 
energy and momentum are not conserved regardless of the timestep with the hybrid appr oach because of the non-conservative source term in )27t . However, 
the solution does not suffer from numerical instabilities (Figure|6] right panel; Rosswog & Price 20071) . 

interested in differences between using the vector potential compared with the Euler potentials, in particular the effect of solving J2 1 1 instead 
of 



4.4 3D Orszag-Tang vortex 

Given that the consistent (energy conserving) vector potential formulation is unstable to the clumping instability for the Orzsag-Tang Vortex 
problem ( ij4.3t — which we expect (and find) is equally true 3D as much as it is in 2D — we consider only the hybrid formulation in 3D, 
by which we mean using equation J27t for the force instead of J42t . The advantage of a 3D problem is that the evolution equation for the 
vector potential, rather than simply being dA z /dt = in 2D, is given by J2U . implemented as d24t > which also differs considerably from the 
evolution equations for the Euler potentials ((2). 
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Figure 8. Results of the Orszag-Tang Vortex evolution in 3D, using 128 X 128 X 16 particles (top, using our NDSPMHD code) 100 X 100 X 10 particles 
(middle, using PHANTOM) and 200 X 200 X 20 particles (bottom, using PHANTOM). Here we have adopted the hybrid vector potential formulation that is 
stable to clumping instabilities and gives good results in 2D. In 3D we observe exponential growth of the A x and A y components of the vector potential 
(right panel, showing the maximum value as a function of time, alongside the evolution of total energy). When these components grow to the same order of 
magnitude as A z , large, low density voids appear in the solution (left panels), together with an exponential divergence in total energy (right panels). 



In 3D we set up the problem similarly to Dolag & Stasvszvr] 1 20081) . by adding a shortened z dimension to the 2D box, placing particles 
initially on a cubic lattice in the domain x, y £ [—0.5, 0.5] and z G [—0.0625, 0.0625]. The setup parameters are identical to the 2D problem 
described in ^4. 3 1 apart from the 3D domain. We have first computed the problem using the same NDSPMHD code that we have used above 
for the 2D problems, using zero artificial resistivity as in the 2D solution shown in lRosswog & Pried ( 120070 . Given our findings below, we 
h ave also computed the solu tion as a consistency chec k using the implementation of the hybrid scheme in our PHANTOM SPH code (used 
in Kitsionas et al. i 20091) and Price & Federrathl (2009)), which being parallelised, was also used to compute a higher resolution version. We 
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show results in Figure[8]using 128 x 128 x 16 particles with NDSPMHD (top row), 100 x 100 x 10 particles using PHANTOM (middle row) 
and 200 x 200 x 20 particles (bottom row), also using PHANTOM. 

The results initially (not shown) are similar to the 2D results — and therefore quite reasonable — but only for a finite time. We find 
that in the 3D case the initially zero A x and A y components of the magnetic field grow exponentially with time, illustrated in the right 
panels of Figure[8] which show the maximum (absolute) value of the A x , A y and A z components of the vector potential as a function of time 
(bottom panel in the right hand figures), alongside the evolution of the total energy (top panel in the right hand figures). The growth of these 
unphysical components is initiated by simple round-off error, evident from the differing starting values between the PHANTOM results (bottom 
two rows) which stores the velocity and magnetic fields as 4-byte variables (giving a starting value of max(A x , A y ) ~ 10 -8 ), compared 
with the NDSPMHD results (top row), which stores both fields to 8-byte precision, giving a starting value of max(A x , A y ) ~ 10" 16 similar to 
the round-off error in (e.g.) momentum conservation (Figure|7]( using double precision variables. The growth rate is exponential in all cases, 
with max(A x , A y ) oc~ e 20t for the top and bottom rows in Figure[8]and max(A x , A y ) oc~ e 10t for the middle row. 

At the point where the A x and A y components become similar in magnitude to the (physical) A z component, large low density voids 
appear in the density field (left panels of Figure [8]l and correspondingly exponential growth in the total energy is observed (right panels 
of Figure [8]l, bringing the simulation rapidly to a halt. We also find the same outcome, though with variations in the exact time at which 
the simulation is disrupted, regardless of whether the Galilean invariant gauge l !22b or the standard gauge {23} is used to evolve the vector 
potential. In general the exact nature of the disruption caused by the growth of unphysical components of the vector potential in 3D depends 
on small details such as round-off error in the code and the nature of the problem studied, however we find similar problems attempting to 
use the vector potential in 3D star formation problems. It should be noted that no such problems arise with the use of the Euler potentials for 
the 3D Orszag-Tang Vortex, since the dA z /dt = evolution corresponds to the da/dt — part of (O whilst the /3 variable is initially set to 
the 2 position of the particles, in which we observe no change, corresponding to the d/3/dt — part of l|2}. 

Actually the results shown in Figure[8]bring our investigation full circle. In fact, we started our examination of the vector potential as an 
alternative to the Euler potentials by studying the 3D Orszag-Tang Vortex, finding the results discussed above. The hope was that a consistent 
formulation of the vector potential from a variational principle, that guarantees exact energy conservation in the code and furthermore 
directly couples both the gauge d — 1 b and the numerical formulation of dA/dt equation d24b to the force formulation d42b would resolve the 
instabilities observed in Figure [8] However, as has already been discussed, the consistent formulation turns out to be itself unstable to the 
tensile instability known to plague conservative formulations of standard SPMHD. The journey is summarised below. 



5 DISCUSSION 

In this paper we have considered the use of the vector potential as a representation for the magnetic field in the context of the Lagrangian 
Smoothed Particle Magnet ohydrodynamics method. In particular, we have addressed the question of whether using the vector potential may 
resolve the issues relating to the restrictions placed on the evolution of 3D magnetic fields using the Euler potentials. To this end we have 
derived a consistent, Hamiltonian formulation for vector potential SPMHD that guarantees both the conservation of momentum and energy. 
The formulation itself relies only on the choice of SPH formulations for the density summation ( i]2.2.U . for obtaining the magnetic field 
from the vector potential via B = V x A ( i)2.2.2b and for the evolution equation for A ( i)2.2.3b . the latter of which involved an appropriate 
choice of gauge, which we required to be Galilean invariant in order to obtain momentum conservation in the equations of motion. 

From these three simple definitions, for which we have used standard variable-smoothing-length SPH operators, we have shown in Q2A 
that the equations of motion can be derived self-consistently from a Lagrangian variational principle, resulting in a form that, reflecting the 
symmetries inherent in the Lagrangian and associated constraint equations, indeed conserves momentum, energy and entropy simultaneously. 
The result is an expression for the MHD (Lorentz) force in SPH that is unique to the vector potential and which differs from all previous 
SPMHD force formulations. The force, given by (142b . or more compactly by d44b can be broken down into components that are non-zero 
either for different numbers of spatial dimensions or depending on whether or not an external magnetic field is applied. 

The expression for the force given by d42b . particularly the "2D" component, initially gave us hope that it might, after all, have been 
possible to formulate equat ions of motion for SPMHD that are both conservative and stable with respect to the tensile or "clumping' 
instability I Monaghan 2000) that prevents the use of exactly conservative force f ormulations in the standard SPMHD approach (B0rv e et al 



2004 ; Paper I ; Paper III ) and was the cause of many initial problems with SPMHD i Phillips & Monaghanll985l) . This was not an unreasonable 



expectation, since the tensile instability in standard SPMHD is caused by non-zero V • B terms in the equations of motion — whereas for the 
vector potential the knowledge that V ■ B = can be "built-in" to the equations of motion which are derived from B = V x A. However 
it turns out that the 2D term is balanced by a highly unstable isotropic term, for which the net pressure is negative (the trigger point for the 
tensile insta bility) when 3/2_B 2 //ip > P, that is, worse than for the standard SPMHD case and not even related solely to the anisotropic part 



of the force l Price [20040 . In numerical tests (Sj4] particularly i]4.2l and i]4.3l > we have found that the instability indeed sets in at around the 



above threshold. Whilst we find that adding correction terms to the force ( i)2.4.5t — though immediately violating the conservation properties 
we so desired — does provide stability for certain problems, we find that the degree of correction required by other problems (i )4.3t starts to 
modify the solution unphysically. In principle a full stability analysis of the consistent vector potential formulation d42b would at least yield 
insights into the exact regime of stability, though we are not confident that much would be gained from such an analysis, let alone a suitable 
correction. 
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The form of the isotropic term causing the tensile instability in d42b arises primarily from the weighti ng with resp ect to p inherent in 
our SPH expression for B = V x A, equation 10, for which we have used a standard curl operator (e.g. Price 2004). It is possible that 
re-formulating this equation using different weightings of p (e.g., using 1/ pi inside the summation) would result in a stable formulation, and 
this would certainly be an approach worth pursuing. However, we also find problems with the vector potential evolution equation even using 
a standard SPMHD force ( i)4.4b . so it is not clear that solving the tensile-instability related problems with the consistent formulation would 
necessarily solve all of the problems. 

An issue not addressed in detail in this paper given the severity of other problems was the fact that the 2D term in J42t involves direct 
second derivatives of the kernel, which are known to be particularly poor using the cubic spline kernel dB2b since it is discontinuous in 
the second derivative. We have side-stepped the issue in this paper by using the smoother quintic kernel dB3b where appropriate, though 
remarkably we find quite reasonable results could be obtained in many cases (for example the ID shock tube problems in !H.l|> even usin g 
the cubic spline. The issue of formulating second derivatives is well known in SPH (see e.g. iBrookshawll 19851 : |Pricel20o4 [Monaghan| [2005h . 
however it is fairly straightforward to show that the usual Brookshaw formulation is equivalent to simply choosing a more appropriate form 
of the kernel for use in the second derivative calculation — that is, using the second derivative of a bell-shaped kernel appropriate to density 
estimates is perhaps not the best approach. For the vector potential we have the freedom to choose the kernel that enters the second derivative 
term, which arises from the kernel in i ] 1 St . completely separately from the kernel used in the density sum. That is, with no loss of consistency, 
the kernels in J 1 3b and d 14b do not have to be the same (we have merely assumed that they were in this paper for simplicity). Thus an 
investigation into the best form of the kernel for computing the curl operation j 1 3b separate from the density sum, perhaps with the second 
derivative estimate J42b in mind, would be a worthwhile exercise. 

Perhaps the most useful aspect of this paper — apart from acting as a warning to the reader intent on similar endeavours — is the 
formulation o f dissipative terms fo r th e vector potential presen ted in [J3] In particular it is clear from this section that the dissipative terms 
formulated bv lPrice & Batd J2007I) and lRosswog & Pried d2007r) for the Euler potentials were not correct, containing V?7 terms which should 
not be present, and more seriously not guaranteeing positive definite dissipation (though both those papers used a term in the thermal energy 
equation that is positive definite, but which instead violates the exact conservation of energy). From an SPH algorithms perspective, Sj3]very 
nicely illustrates the conjugate relationship between the standard J 1 8b and symmetric ( 120b SPH curl operatorQ Using the standard curl for 
A necessitates the use of the symmetric curl for J in order to obtain both energy conservation and positive definite dissipation. We have also 
shown in fJ3]how the resistivity parameter r\ may be constructed appropriate to an artificial resistivity term, demonstrating that this approach 
works well on standard shock tube problems ( §4. lb where dissipation is important. 

In terms of finding a satisfactory approach to maintaining the divergence constraint in three dimensional SPMHD simulations without 
the restrictions associated with the Euler Potentials formulation, we intend to examine generalised forms of the Euler Potentials which can 
represent arbitrary MHD fields and which can be more easily adapted for non-ideal MHD. This is deferred to a future paper. 



6 CONCLUSION 

In summary, we find that using a hybrid formulation of the vector potential in SPH, evolving A using d24b and calculating the force using 
one of the standard stable SPMHD force expressions d27b or {25}, is the approach with the fewest difficulties, and gives good results on one 
and two dimensional test problems - identically to the Euler potentials in 2D. However, we also find problems with this approach for 3D 
problems, leading us to conclude that use of the vector potential is not a viable approach for SPMHD. 
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APPENDIX A: PERTURBATION TO THE MAGNETIC FIELD 

The derivation of the Lagrangian perturbation for the magnetic field, discussed in ^2.4. H and used in equation d32b . is easier to understand if 
one considers that the Lagrangian perturbation S = A + <5x • V is similar to taking a Lagrangian time derivative d/dt = d/dt + v ■ V (more 



4 This has been noted previously bv lCummins & Rudmanl (l999). Here we have demonstrated the equivalent relationship for the variable smoothing length 
SPH operators. 
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precisely d/dt = 8 /St). Taking the time derivative of d!8t . assuming dB ex t/dt = 0, gives 

= — i— m b (A„ - Ab) X ^ [Vq.Wa.fr (ftq)] 

Bi„t rfPa Bint df2 a 



+ 



^E 



dt 



dt 



pa dt Q a dt 



The time derivative of the kernel gradient can be shown (Appendix|Bj to be given by 

d V7W7TJ/ fU \ , dVWab(h a ) dha 

Tt (y«w ab ) = (v ab . v)vw ab (h a ) + dha — , 

the latter term arising only in the case of a variable smoothing length. Using this expression in dAlt and assuming h — h(p) gives 

^f- = T^— y2m b (A a -A b ) X [(v a -Vfr)- V]V tt Wa&(Aa) 

D 

1 v— v /dA a dA 6 \ B int dp a B int dQ, a H a dp a 

+ ^E"M — - — ) xv.^(M-- ¥ "—— + —— < 



(Al) 



(A2) 



Pa di fia p a dt 



where we have defined 

dh a 1 

<9pa fia 



a/la 1 IK A \ 9V a Wai,(/l a ) 
H a EE — } m b A a - A 6 X — — * '-. 



(A3) 



(A4) 



Replacing v with 8~x./5t and time derivatives d/dt with 5 /St we have 
5B a = — y^m b (A a — A b ) x [(<Sx a — 5x. b ) ■ V]V a W ab (h a ) 

"aPa ^ 

— m b (SA a - 5A b ) X VaWab(ha) - ^^-8p a - ^P^5Q a + — <5pa, 



OaPa 



Pa 



Pa 



This is an SPH expression for the Lagrangian perturbation of the magnetic field, i.e., 
SB = V x [<5A - (5x ■ V)A] + <5x • V(V x A). 

The perturbation required in £12.4. 1 1 i.e., 5(p b B b ) = pb<5B + 5p b B b is thus given by 

<5(p D Bi,) = — m c (A 6 - A c ) x [(<5x - 5x e ) ■ V] X7 b W bc (h b ) 

c 

+ ^j- ^ m c (<5A 6 - 5A C ) x V 6 Wb c (frb) + B ea;t <$p6 - B ''-"' t P b jq 6 + H 6 <5p ( 



(A5) 



(A6) 



(A7) 



The higher order perturbations <5fi relating to the smoothing length gradients are of order the truncation error of the SPH method (i.e., 
0(h 2 )) and thus may be justifiably neglected from the analysis. On the other hand, to obtain a method that conserves energy to timestepping 
accuracy it is necessary to include these terms. Thus, for completeness, the perturbation to f2, from Equation d!9t , is given by 

(a , dh a u X X \ V71 d W ab {ha) 



SQ a = — — Sp a - — - m b \(5xa - <5x b ) • V 
p a a p a ^ 

where we have defined the dimensionless quantity 



dha 



Co. = Pa 



8 2 hg 

dpi 



E 



in b 



dWg b (h g) , (dh. 

dh a 



+ E 



m b 



d 2 Wg b (h a 
dh\ 



(A8) 



(A9) 



6 " v-f"/ b 

Combining dA8t and dA7b . the full perturbation including all terms relating to the smoothing length gradients is given by 



5(p b B b ) = -j- ^ rn c (A b - A c ) x [(<5x D - 5x c ) ■ V] X7 b W bc (h b 

c 

+ 7^- E mc ( SAb _ ^ Ac ) X ^bW bc {h b ) + Bext5p b 



+ 



H. B int > 
b H f; — C,t> 



B b , int p b dh b t-^ dW bc {h b ) 
0Pb H p; TT— 2^ "ic [(ox 6 - <5x c ) ■ V D ] 



Hb <9pb 



5/ib 



(A10) 



APPENDIX B: KERNEL FUNCTIONS AND DERIVATIVES 

In this appendix we describe the various kernel functions used in this paper and derive the derivatives required in the vector potential 
formulation, showing how they may be computed from the dimensionless kernel functions. 
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Bl Kernel function 



It is common to express the SPH kernel in the form 

W ab (h) EE 



(Bl) 



where a is a normalisation constant, h is the smoothing length, v is the numbe r of spatial dimen sions and /(g) is a dimensionless function 
of the variable q = \r a — *b\/h. For the standard SPH cubic spline kernel (e.g. Mon aghar]| 1992 ) the function / is given by 



/(<?) = 



±(2-g) 3 -(l-g) 3 , 0<«<1; 
|(2 -g) 3 , l^g<2; 
0, g ^ 2. 



(B2) 



with normalisation <j = [2/3, 10/ (7ir), 1/vr] in [1, 2, 3] dimensions. In this paper we also use the quintic spline kernel given by (e.g. Morris! 
1 1996b 



/(«) = 



'(3-g) 5 -6(2-g) 5 + 15(l-g) 5 , < g < 1; 

(3-g) 5 -6(2-g) 5 , 1 < g < 2; 

(3-g) 5 , 2s;g<3; 

0, g > 3. 



(B3) 



with normalisation a — [1/24, 96/(11997r), 1/(20tt)]. The advantage of the quintic is that it more closely approximates the Gaussian by 
extending the compact support radius to 3h and has smooth second derivatives — important here since we use the second derivative directly 
in the force (equationl42b. The disadvantage is that it is considerably more expensive to compute, particularly in three dimensions, where the 
neighbour number, and thus the cost, increases by a factor of (3/2) 3 sa 3.4. 



B2 Kernel first derivatives 

The derivative of the kernel with respect to the particle coordinates, holding the smoothing length constant is given by 

V a W ab = £f{q)Vq = tabj^f'iq), (B4) 

where we have used the fact that Vg = V(|r a (,|//i) = r a b/h, defining r a t = r a — r b . The kernel derivative can also be written in the form 



VaWab = tabFab, where F ab ee j—^f'(q). 



(B5) 



Note that the definition we use for F ab differs slightly from that used by iMonaghanl < 19921) since we use the unit vector in our definition, 
reflecting the implementation in our code. 

The derivative of the kernel with respect to h is given in terms of /(g) by 

dWab o ,, s , a i dq a r , -, 



(B6) 



B3 Second derivatives 

The second derivative of the kernel with respect to particle coordinates is given by 

viviwab = vi [fi^fiq)] , 

1 s>, 



h u+2 

For the Laplacian this reduces to 



V'Wab 



h v+2 



f"(q) + (v-l)±f'(q) 



The time derivative of the kernel gradient can be derived from equation JB4b . giving 

d , Vabjrab ■ VWqi,) (v aft ■ X7W ab )r ab (V ab ■ V ab ) af"(q) _ r^yyu, 

di iVaWab) = w~bj 2 m 1 ab = ( ab ' } ab ' 



The mixed second derivative of the kernel with respect to particle coordinates and h is given by 



-^-(VaWab) = -iy + i)^+2f'{q)iab + j^f"{q)^T a b, 



F„ h r„ 



qf"(q)Kb- 



(B7) 
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Finally, the second derivative with respect to the smoothing length is given by 

d 2w ab i I i \ cr ./ x . cr g , erg 2 „ 

-^r = + 1)^5/(9) + + + 2)1 if W + j^f («). 

_ u(v + l) 2Q+l)g aq 2 „ 

- ^2 ^ ^ ab + J^J W- (BH) 



APPENDIX C: TRANSLATION OF THE SPH EQUATIONS OF MOTION IN THE CONTINUUM LIMIT 

The MHD equations of motion are given by 
dv VP J x B 



(CI) 



dt p p 

In this appendix we show that the SPH equations of motion derived in [J2A4] translate to the above in the continuum limit. Written in terms 
of the vector potential, ( IClt becomes 



dv _ VP J x (V x A) J x B ext 
dt p p p 

or, in tensor notation, 



(C2) 



dv_ 

dt 



1 dP + 1 

p dx 1 p 



J- 



dA 3 



J 



OA' 



P R fc 



(C3) 



dx 1 dx 3 

The terms in the SPH equations of motion d42l >, d44b or d47| > can be decoded into continuum form using the SPH summation interpolant 



^ m b — Wab A a , 



(C4) 



the derivatives of which give 

A b dWab dA a 



E 



m b 



pb dxi 



UI± a \ - 

~dxi'' ^ 



A b d 2 W ab d 2 A a 



pb dxixi dxixl 



(C5) 



In translating the equations of motion we firstly neglect all terms relating to gradients in the smoothing length - that is, we provide the 
translation of l !47t rather than (142b . Whilst in principle the extra terms in J42b could also be translated, the proof that they are correctly 
derived lies in the fact that the numerical equations demonstrate the conservation of energy to timestepping accuracy, which we have shown 
in £0 For a constant smoothing length, the SPH equations of motion J47t can be written in tensor notation (akin to equation l44l in the form 



dt 



E 



nit 



Pi 



si 

Pi 



dW ab 

dxi 



po V 



A\) 



Bi + Bi 

Pa Pi 



dxidx l n 



(C6) 



where 



-P8 13 + 



po 



B l Bi xt 



+ 5 l - 



1°' 



2B B, 



A 1 J 3 



(C7) 



The first term in l lC6b is similar to the usual conservative SPMHD force (e.g. Paper II), albeit with a different stress tensor, and can be 
straightforwardly translated to 1/p dS 13 /dx 3 . Expanding this using dC7b . we have 



~~dt 



— + B 

dx % 



dx 3 dx 1 



2PLt 



dB 3 rl dA r Al dJ 3 
J — — - — A 



+ 2nd term 



dx 1 dx 3 dx 3 
The second or "2D term" in dC6l > can be translated using dC5b as follows 

'•' - ^ r J^_ _ €juB 3 ' d 2 (pA k ) <-,,< i> 

po " " ' ' 



2nd term = — - — A 



po 



p 2 dx i dx 



' d 




dx 1 





+ 



d 

dx 1 



A B 3 



(C8) 
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po P 2 dx i dx l ' po dx % 

where in the third term we obtain a change of sign because the indices j and k have changed order with respect to the Levi-Civita tensor. 
Expanding the derivative terms and collecting/cancelling terms where appropriate, after some straightforward algebra we are left with 

k .3 T33 a a k r M f 

(CIO) 



e jfei B 3 d 2 A k e jkl 8A k dB 3 e m dA k dB 3 
* « ; r H ^ r o „■ r 



po p dx i x l pop dx 1 dx 1 



pop dx i dx 1 



which, using B 3 int = e ]k idA k /dx 1 and J k = 1/ p e k jidB 3 /dx 1 gives 



o B 3 0B 3 nt 
pop dx 1 



B 3 nt dB 3 | 18A k jk 
pop dx 1 p dx" 1 



(Cll) 
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Putting this together with JC8b and collecting terms, noting that _B J = B\ nt + B 3 ext and that we have previously assumed dB 3 , xt /dx 1 ' — 
we have 



~~dt 



dP 

dx l 



Mo 



dB i _ dB J 
dxi dx l 



dA J 
dx l 



J 3 — J 



dA l 
~dxi 



,idJ j 
i 

dxi> 



(C12) 



The last term is zero in the continuum limit since it is the divergence of a curl (V • J) = V • (V X B)//ito. This completes the proof, giving 
the equations of motion in the form 



~dt 



dP [ Bj xt 

dx i fAQ 



dB l dB j 
dxi Ox 1 



dx l dxi 



(C13) 



which, upon expanding [(V x B) x B ext ] z = B J ext (dB i /dx j - dB j /dx 1 ), is identical to < fC3T >. 
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